library(ggplot2)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
setwd("/Users/adriancerezuelahernandez/Desktop/Q4/AD/SerTemp/Projecte")
ser=ts(read.table("/Users/adriancerezuelahernandez/Desktop/Q4/AD/SerTemp/Repositori de series temporals/ConsumElec.dat"),start=1990,freq=12)
autoplot.zoo(as.zoo(ser), facets=NULL) + aes(colour = NULL, linetype = NULL) + facet_free() + ggtitle("Consumo de electricidad(energía final) en España") + xlab("Tiempo") + ylab("GwH") + geom_point()
m <- apply(matrix(ser,nrow=12),2,mean)
## Warning in matrix(ser, nrow = 12): data length [347] is not a sub-multiple or
## multiple of the number of rows [12]
v <- apply(matrix(ser,nrow=12),2,var)
## Warning in matrix(ser, nrow = 12): data length [347] is not a sub-multiple or
## multiple of the number of rows [12]
plot(v~m)
boxplot(ser~floor(time(ser)))
La varianza no es constante, será necesario aplicar logaritmos a la
serie.
lnser = log(ser)
autoplot.zoo(as.zoo(lnser), facets=NULL) + aes(colour = NULL, linetype = NULL) + facet_free() + ggtitle("Consumo de electricidad(energía final) en España") + xlab("Tiempo") + ylab("GwH(log)") + geom_point()
Realizando de nuevo el análisis de la varianza para la serie
transformada
m <- apply(matrix(lnser,nrow=12),2,mean)
## Warning in matrix(lnser, nrow = 12): data length [347] is not a sub-multiple or
## multiple of the number of rows [12]
v <- apply(matrix(lnser,nrow=12),2,var)
## Warning in matrix(lnser, nrow = 12): data length [347] is not a sub-multiple or
## multiple of the number of rows [12]
plot(v~m)
boxplot(lnser~floor(time(lnser)))
Ahora se ve que se ha homogenizado la varianza, y es constante. En el
boxplot se puede observar también la presencia de algunas observaciones
atípicas, las cuales serán tratadas adecuadamente más tarde.
Análisi del patrón estacional
monthplot(lnser)
ts.plot(matrix(lnser,nrow=12),type="o",col=1:8)
## Warning in matrix(lnser, nrow = 12): data length [347] is not a sub-multiple or
## multiple of the number of rows [12]
Subidas de consumo en meses centrales como Enero y Julio, la serie
presenta un patrón estacional. Será necesario aplicar una diferenciación
estacional, y trabajar con un incremento de la estación o mes respecto
del año anterior.
d12lnser = diff(lnser,lag=12)
plot(d12lnser)
abline(h=0,col=2)
Claramente, a través del gráfico de la serie se puede ver como la media
no es constante. Por tanto, será necesario aplicar tantas
diferenciaciones regulares como sea necesario con tal de corregirla. Una
vez se pueda considerar que es prácticamente constante, deberá verse qué
diferenciación es la que da una mejor relación entre mantener la media
constante y no aumentar la varianza en exceso.
d1d12lnser = diff(d12lnser)
plot(d1d12lnser)
abline(h=0,col=2)
d1d1d12lnser = diff(d1d12lnser)
plot(d1d1d12lnser)
abline(h=0,col=2)
var(lnser)
## V1
## V1 0.05764872
var(d12lnser)
## V1
## V1 0.00200361
var(d1d12lnser)
## V1
## V1 0.00219549
var(d1d1d12lnser)
## V1
## V1 0.006041761
La varianza de la serie tras la última diferenciación regular ha aumentado considerablemente. Tras aplicar una sola diferenciación regular también aumenta la varianza, pero es un aumento mínimo, por lo que se puede tomar a costa de hacer constante la media. Por tanto, la última diferenciación regular no se tendrá en cuenta y, tras realizar las transformaciones, se trabajará con la serie estacionaria \(W_t\), que es de la forma \[W_t = (1-B)(1-B^{12})log(X_t)\]
A continuación, se presentará el ACF (Auto-Correlation Function) y el PACF(Partial Auto-Correlation Function) de la serie transformada, con el fin de decidir cuáles serán los parámetros de los modelos a presentar.
acf(ser,ylim=c(-1,1),lag.max = 72) #Para ver que la serie original no es estacionaria
par(mfrow=c(1,2))
acf(d1d12lnser, ylim=c(-1,1), lwd=2, lag.max=72,col=c(2,rep(1,11)))
pacf(d1d12lnser, ylim=c(-1,1), lwd=2, lag.max=72, col=c(rep(1,11),2))
par(mfrow=c(1,1))
Parte estacional: MA(1) estacional. Para el PACF se observan dos retardos significativos. Después del tercero tenemos dos significativos, pero como sobresalen por poco de las bandas de confianza, se considerarán valores aleatorios fruto del azar, y el modelo que tomaremos será un AR(2) estacional.
Parte regular: En el ACF encontramos tres retardos significativos. Los posteriores, como antes, se consideraran fruto del azar, considerando también que están muy lejos del inicio. Se considerará un MA(3) regular en este caso. Para el PACF se observan cinco u ocho retardos significativos, según si consideramos los últimos dos valores aleatorios. En un principio se propondrá un AR(8) regular, con la opción de poder reducir los parámetros más tarde en la estimación.
Combinando las opciones consideradas, y teniendo en cuenta que previamente se ha aplicado una diferenciación regular y una estacional, los modelos iniciales propuestos serían los siguientes: \[ARIMA(0,1,3)(0,1,1)_{12}\] \[ARIMA(0,1,3)(2,1,0)_{12}\] \[ARIMA(8,1,0)(0,1,1)_{12}\] \[ARIMA(8,1,0)(2,1,0)_{12}\]
De cara a la estimación de modelos se tomarán los dos modelos con un MA(1) estacional.
(mod1 <- arima(d1d12lnser, order=c(0,0,3), seasonal=list(order=c(0,0,1),period=12)))
##
## Call:
## arima(x = d1d12lnser, order = c(0, 0, 3), seasonal = list(order = c(0, 0, 1),
## period = 12))
##
## Coefficients:
## ma1 ma2 ma3 sma1 intercept
## -0.5682 -0.1161 -0.0774 -0.8006 -1e-04
## s.e. 0.0557 0.0633 0.0578 0.0390 1e-04
##
## sigma^2 estimated as 0.0009004: log likelihood = 690.59, aic = -1369.19
Es preferible que la media estimada de la serie transformada no sea significativamente diferente de 0, indicador de que la serie transformada está centrada alrededor de 0. Según el ratio del intercept dado por el modelo ajustado, es el caso. De este modo, se toma la serie sin diferenciaciones, ni regular ni estacional, y se ajustan con el modelo a continuación.
(mod1 <- arima(lnser, order=c(0,1,3), seasonal=list(order=c(0,1,1),period=12)))
##
## Call:
## arima(x = lnser, order = c(0, 1, 3), seasonal = list(order = c(0, 1, 1), period = 12))
##
## Coefficients:
## ma1 ma2 ma3 sma1
## -0.5672 -0.1152 -0.0758 -0.7942
## s.e. 0.0559 0.0632 0.0577 0.0383
##
## sigma^2 estimated as 0.0009034: log likelihood = 690.22, aic = -1370.43
Si hacemos el test de t-ratios para cada coeficiente se podrá comprobar si es conveniente rechazar alguno y establecerlo a 0 en el modelo.
cat("\nT-ratios", round(mod1$coef/sqrt(diag(mod1$var.coef)),2))
##
## T-ratios -10.15 -1.82 -1.31 -20.74
Según el test realizado, los coeficientes \(ma2\) y \(ma3\), tienen un valor absoluto menor que dos. Por tanto, deben fijarse a 0 en el modelo a continuación. No obstante, \(ma2\) es dudoso, pues la diferencia con el valor umbral es muy poca, por lo que puede ser preferible mantenerlo en el modelo. De hecho, si comprobáramos el AIC del modelo fijando este parámetro a 0, el resultado sería un modelo peor que el planteado a continuación.
(mod1 <- arima(lnser, order=c(0,1,3), seasonal=list(order=c(0,1,1),period=12),fixed=c(NA,NA,0,NA)))
##
## Call:
## arima(x = lnser, order = c(0, 1, 3), seasonal = list(order = c(0, 1, 1), period = 12),
## fixed = c(NA, NA, 0, NA))
##
## Coefficients:
## ma1 ma2 ma3 sma1
## -0.5861 -0.1636 0 -0.7924
## s.e. 0.0522 0.0489 0 0.0383
##
## sigma^2 estimated as 0.0009085: log likelihood = 689.36, aic = -1370.72
Se puede ver como la mejora después de fijar los dos parámetros es muy ligera, prácticamente mínima, respecto al modelo completo.
cat("\nT-ratios", round(mod1$coef/sqrt(diag(mod1$var.coef)),2))
## Warning in mod1$coef/sqrt(diag(mod1$var.coef)): longer object length is not a
## multiple of shorter object length
##
## T-ratios -11.23 -3.35 0 -15.18
Al volver a calcular los ratios de cada coeficiente, para comprobar si se ven afectados, se comprueba que todos siguen siendo significativos, por lo que no se podrá fijar a 0 ningún otro.
Una vez estimado el primer modelo, a continuación se estima el segundo, el cual contiene un AR(8) regular en este caso.
(mod2 <- arima(d1d12lnser, order=c(8,0,0), seasonal=list(order=c(0,0,1),period=12)))
##
## Call:
## arima(x = d1d12lnser, order = c(8, 0, 0), seasonal = list(order = c(0, 0, 1),
## period = 12))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.5792 -0.4566 -0.4525 -0.3917 -0.2753 -0.1726 -0.1684 -0.1125
## s.e. 0.0554 0.0633 0.0667 0.0693 0.0697 0.0669 0.0633 0.0554
## sma1 intercept
## -0.7964 -1e-04
## s.e. 0.0400 1e-04
##
## sigma^2 estimated as 0.0008851: log likelihood = 693.48, aic = -1364.96
Igual que para el modelo anterior, el ratio del intercept es menor que |2|, por lo que la media no es significativamente diferente, y se toma la serie sin diferenciaciones previas realizadas, ajustándolas con el propio modelo.
(mod2 <- arima(lnser, order=c(8,1,0), seasonal=list(order=c(0,1,1),period=12)))
##
## Call:
## arima(x = lnser, order = c(8, 1, 0), seasonal = list(order = c(0, 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.5787 -0.4554 -0.4506 -0.3889 -0.2732 -0.1709 -0.1676 -0.1116
## s.e. 0.0554 0.0634 0.0668 0.0693 0.0698 0.0670 0.0635 0.0554
## sma1
## -0.7913
## s.e. 0.0394
##
## sigma^2 estimated as 0.0008874: log likelihood = 693.19, aic = -1366.37
Si hacemos el test de t-ratios para cada coeficiente se podrá comprobar si es conveniente rechazar alguno y establecerlo a 0 en el modelo.
cat("\nT-ratios", round(mod2$coef/sqrt(diag(mod2$var.coef)),2))
##
## T-ratios -10.44 -7.19 -6.75 -5.61 -3.91 -2.55 -2.64 -2.01 -20.07
Todos ellos están por encima de 2, en valor absoluto. Por lo tanto, ningún coeficiente del modelo deberá ser fijado a 0.
En primer lugar, se comprobará a través del ACF y PACF de los residuos, que todas las observaciones caen dentro de las bandas de confianza, tal como se espera que lo hagan.
resi=resid(mod1)
par(mfrow=c(1,2))
acf(resi,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
pacf(resi,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
par(mfrow=c(1,1))
plot(resi)
abline(h=0)
abline(h=c(-3*sd(resi),3*sd(resi)),lty=3,col=c(1))
Tal como se esperaba, todas las observaciones, a excepción de un par en
el PACF, caen dentro de las bandas de confianza. El mencionado retardo
que se sale de estas, lo hace por muy poco y podría considerarse un
valor aleatorio fruto del azar.
Al realizar una primera exploración de los residuos del primer modelo, se observa que, sin tener en cuenta el claro pico pronunciado, y un segundo pico que se sale de las bandas de confianza, que aparecen entre los años 1990 y 1995, la variabilidad es similar para toda la serie. Por lo tanto, se puede concluir que este no es un caso de varianza no constante producida por una alta volatilidad en los datos, sino de la presencia de observaciones atípicas.
De todas maneras, otra representación gráfica que puede ayudar a decidir frente a qué situación nos encontramos es el scatterplot de la raíz cuadrada del valor absoluto de los residuos.
scatter.smooth(sqrt(abs(resi)),lpars=list(col=2))
Vemos que, globalmente, la línea roja es horizontal y los residuos se
distribuyen de manera homogénea alrededor de ella, a excepción de
algunas observaciones atípicas, que tal como se había visto en el
anterior gráfico, se sitúan entre 1990 y 1995. Por lo tanto, la
conclusión es la misma, la variancia de los residuos es constante.
A continuación se estudiará su normalidad. En primer lugar, el gráfico de normalidad de los residuos es el siguiente:
qqnorm(resi)
qqline(resi,col=2,lwd=2)
Se observa como en los extremos aparecen colas que se alejan de la
recta, pero no lo suficiente como para estar hablando de volatilidad.
Son debidas a la ya mencionada presencia de atípicos. Se puede suponer,
de momento, que se cumple la hipótesi de normalidad.
Mediante el histograma con la curva normal superpuesta y el test de Shapiro-Wilks se comprobará esta suposición.
hist(resi,breaks=20,freq=FALSE)
curve(dnorm(x,mean=mean(resi),sd=sd(resi)),col=2,add=T)
shapiro.test(resi)
##
## Shapiro-Wilk normality test
##
## data: resi
## W = 0.98814, p-value = 0.006275
Todo y ser el histograma únicamente un gráfico ilustrativo, se puede ver que a priori se debería de cumplir la hipótesis de normalidad. En cuanto al test de Shaprio-Wilks, el p-value tiene un valor de 0.0006275, que es inferior a 0.05 y por tanto rechaza la hipotesis nula, concluyendo así que no se puede garantizar que se cumpla la hipótesis de normalidad para los residuos del modelo.
Para finalizar el análisis, se estudiará la independencia de estos y si se trata o no de ruido blanco, a través del test de Ljung-Box.
tsdiag(mod1,gof.lag=72)
Para que se pueda hablar de independencia y ruido blanco, todas las
observaciones en el último gráfico deberían quedar por encima de la
línea discontinua de color azul. Se puede ver como es el caso, por tanto
no se rechaza la hipótesis y concluimos que el modelo \(mod1\) explica suficientemente bien
nuestros datos.
A continuación, se seguirá el mismo procedimiento para el segundo modelo propuesto, empezando, de nuevo, por comprobar que las observaciones en el ACF y PACF caigan dentro de las bandas de confianza
resi=resid(mod2)
par(mfrow=c(1,2))
acf(resi,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
pacf(resi,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
par(mfrow=c(1,1))
plot(resi)
abline(h=0)
abline(h=c(-3*sd(resi),3*sd(resi)),lty=3,col=c(1))
Se puede ver como tanto el ACF y PACF, como el gráfico que permite
realizar una primera exploración de los residuos, son muy similares a
los mostrados para el primer modelo, por lo cual se pueden sacar las
mismas conclusiones: no es un caso de varianza no constante producidad
por una alta volatilidad en los datos, sino de la presencia de
observaciones atípicas.
scatter.smooth(sqrt(abs(resi)),lpars=list(col=2))
El mismo caso tenemos con el scatterplot de la raíz cuadrada del valor
absoluto de los residuos, es muy similar al anterior. Por lo tanto, la
conclusión es la misma, la variancia de los residuos es constante.
qqnorm(resi)
qqline(resi,col=2,lwd=2)
hist(resi,breaks=20,freq=FALSE)
curve(dnorm(x,mean=mean(resi),sd=sd(resi)),col=2,add=T)
shapiro.test(resi)
##
## Shapiro-Wilk normality test
##
## data: resi
## W = 0.98639, p-value = 0.002389
En el caso del estudio de la normalidad, en este caso tampoco se puede garantizar que se cumpla. El gráfico de normalidad de los residuos muestra como forman esa línea recta esperada, a excepción de las colas en los extremos por la presencia de atípicos. Mientras tanto, una vez más, el test de Shapiro-Wilks tiene un p-value por debajo de 0.05, de modo que se rechaza la hipótesis nula de que los residuos cumplan normalidad.
Para finalizar el análisis, se estudiará la independencia de estos y si se trata o no de ruido blanco, a través del test de Ljung-Box.
tsdiag(mod2,gof.lag=72)
Igual que anteriormente, se puede hablar de ruido blanco e independencia
en las observaciones. En este caso, tal como se puede observar en el
tercer gráfico, todas están muy alejadas de la línea discontinua, hecho
que prueba que los resultados de este test son aún mejores que los del
primer modelo. Por tanto, se puede concluir que \(mod2\) también es un modelo que explica
suficientemente bien nuestros datos.
A continuación se estudiará la causalidad e invertibilidad del modelo, pero antes se debe tener en cuenta que se decidirán dichas cualidades en función de si se trabaja con la parte MA o AR del modelo. Se dará por hecho que no hay ningún problema si todos los puntos (raíces del polinomio) caen dentro del circulo unidad.
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
plot(mod1)
Se observa que todos los puntos caen en el centro del círculo unidad, así que no vemos ningún indicio de que haya ningun problema.
Se prosigue el estudio con el análisis de la invertibilidad. Se necesitarán únicamente las raíces complejas del polinomio resultante, si los módulos de dichas raíces son mayores que 1 (>1) se podrá garantizar la invertibilidad.
polyroot(c(1,mod1$model$theta))
## [1] 0.5097910+0.8829839i -0.8829839+0.5097910i -0.5097910-0.8829839i
## [4] 0.8829839-0.5097910i 0.0000000+1.0195820i -1.0195820+0.0000000i
## [7] 0.0000000-1.0195820i 1.0195820+0.0000000i -0.5097910+0.8829839i
## [10] -0.8829839-0.5097910i 0.5097910-0.8829839i 0.8829839+0.5097910i
## [13] 1.2618033-0.0000000i -4.8446723-0.0000000i
print("Raíces del polinomio característico")
## [1] "Raíces del polinomio característico"
Mod(polyroot(c(1,mod1$model$theta)))
## [1] 1.019582 1.019582 1.019582 1.019582 1.019582 1.019582 1.019582 1.019582
## [9] 1.019582 1.019582 1.019582 1.019582 1.261803 4.844672
Se concluye que se puede garantizar que el modelo es invertible, debido a que se cumple el requisito anteriormente mencionado.
Para el estudio de la causalidad, los módulos de las raíces deberán de, nuevamente, ser mayores que 1 (>1) para afirmar que el modelo en cuestión sea causal.
polyroot(c(1,mod1$model$phi))
## complex(0)
print("Raíces del polinomio característico")
## [1] "Raíces del polinomio característico"
Mod(polyroot(c(1,-mod1$model$phi)))
## numeric(0)
Dada la salida, se puede ver como no tiene polinomio característico, al tratarse de un modelo MA, tanto estacional como regular. Por tanto, se concluye que el modelo es también causal.
Seguidamente se procede a encontrar la expresión del modelo como AR y MA infinitos. Los coeficientes devueltos por la función \(ARMAtoMA()\) representan sus términos, y describen cómo los errores de predicción anteriores se utilizan para hacer predicciones futuras. Cada coeficiente representa el efecto de un error de predicción específico en la predicción actual.
##Pesos Pi (AR infinit)
(pis=-ARMAtoMA(ar=-mod1$model$theta,ma=-mod1$model$phi,lag.max=36))
## [1] -0.58610426 -0.50710339 -0.39309343 -0.31334833 -0.24795905 -0.19658900
## [7] -0.15578408 -0.12346476 -0.09784719 -0.07754566 -0.06145619 -0.84108646
## [13] -0.50301768 -0.43241009 -0.33572363 -0.26750493 -0.21170519 -0.16784116
## [19] -0.13300425 -0.10541068 -0.08353918 -0.06620629 -0.05246956 -0.66945130
## [25] -0.40095150 -0.34451169 -0.26750950 -0.21314546 -0.16868605 -0.13373505
## [31] -0.10597722 -0.08399077 -0.06656365 -0.05275289 -0.04180752 -0.53064432
La salida muestra los coeficientes del modelo expresado como un AR infinito. Para determinar si es estacionario y/o invertible, deben examinarse estos. Será estacionario si y solo si la suma de los coeficientes al cuadrado converge a un valor finito. Si la suma diverge, entonces no será estacionario. De igual manera se comprueba la invertibilidad del mismo, solo que la suma en este caso es de los coeficientes en valor absoluto. De igual manera, será invertible si esta converge. Si el resultado fuera infinito, en cualquiera de los dos casos, sería necesario tomar medidas con tal de arreglar este problema.
sum(pis^2)
## [1] 3.659775
sum(abs(pis))
## [1] 9.039941
Las dos sumas convergen, se puede afirmar pues que el modelo expresado como AR infinito es estacionario e invertible
##Pesos Psi (MA infinit)
(psis=ARMAtoMA(ar=mod1$model$phi,ma=mod1$model$theta,lag.max=36))
## [1] -0.5861043 -0.1635852 0.0000000 0.0000000 0.0000000 0.0000000
## [7] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 -0.7923814
## [13] 0.4644181 0.1296219 0.0000000 0.0000000 0.0000000 0.0000000
## [19] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [25] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## [31] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
Procediendo de igual manera para un MA infinito, los coeficientes del modelo expresado como tal son los que se muestran a la salida. Se comprueban, también de igual manera, la estacionariedad e invertibilidad del mismo.
sum(psis^2)
## [1] 1.230633
sum(abs(psis))
## [1] 2.136111
Las dos sumas convergen a un valor finito. Por tanto, se puede concluir que el modelo será estacionario e invertible, tanto expresado como un AR infinito, como expresado como un MA infinito.
En cuanto a las medidas de adecuación, se analizan los criterios AIC y BIC, que explican cuán bueno es el ajuste del modelo.
print(c("AIC del modelo 1:",mod1$aic))
## [1] "AIC del modelo 1:" "-1370.7152558585"
print(c("BIC del modelo 1:",BIC(mod1)))
## [1] "BIC del modelo 1:" "-1355.47069188659"
De la misma manera que se ha procedido con el primer modelo, \(mod1\), se estudiará también la invertiblidad, causalidad, expresión como AR y MA infinitos, y medidas de adecuación del segundo modelo, \(mod2\).
Se comenzará, en primer lugar, con la comprobación de que todas las raíces del modelo caigan dentro del círculo unidad.
plot(mod2)
Vemos nuevamente que todos los puntos caen en el centro del círculo unidad, así que no vemos ningún indicio de que haya ningun problema.
Se prosigue el estudio con el análisis de la invertibilidad y causalidad.
polyroot(c(1,mod2$model$theta))
## [1] 0.5098473+0.8830814i -0.8830814+0.5098473i -0.5098473-0.8830814i
## [4] 0.8830814-0.5098473i 0.0000000+1.0196946i -1.0196946-0.0000000i
## [7] 0.0000000-1.0196946i 1.0196946+0.0000000i -0.5098473+0.8830814i
## [10] -0.8830814-0.5098473i 0.5098473-0.8830814i 0.8830814+0.5098473i
print("Raíces del polinomio característico")
## [1] "Raíces del polinomio característico"
Mod(polyroot(c(1,mod2$model$theta)))
## [1] 1.019695 1.019695 1.019695 1.019695 1.019695 1.019695 1.019695 1.019695
## [9] 1.019695 1.019695 1.019695 1.019695
polyroot(c(1,mod2$model$phi))
## [1] 0.7118067-0.000000i -1.1249832+0.973337i -0.0818620-1.361681i
## [4] 0.8873392-1.074464i -0.0818620+1.361681i -1.5737250-0.000000i
## [7] -1.1249832-0.973337i 0.8873392+1.074464i
print("Raíces del polinomio característico")
## [1] "Raíces del polinomio característico"
Mod(polyroot(c(1,-mod2$model$phi)))
## [1] 1.265513 1.317305 1.317305 1.265513 1.283098 1.399223 1.283098 1.399223
Se observa que, de la misma manera que para el primer modelo, todos los módulos son de valor superior a 1, tanto en el caso de las phis como en el de las thetas, por lo tanto se concluye que el segundo modelo es invertible y causal.
Estudio de su estacionariedad e invertibilidad. La expresión de este segundo modelo como AR infinito tiene por coeficientes los siguientes:
##Pesos Pi (AR infinit)
(pis=-ARMAtoMA(ar=-mod2$model$theta,ma=-mod2$model$phi,lag.max=36))
## [1] -0.57865512 -0.45541621 -0.45057650 -0.38890235 -0.27321652 -0.17092807
## [7] -0.16755656 -0.11163512 0.00000000 0.00000000 0.00000000 -0.79133147
## [13] -0.45790800 -0.36038518 -0.35655537 -0.30775067 -0.21620483 -0.13526076
## [19] -0.13259278 -0.08834039 0.00000000 0.00000000 0.00000000 -0.62620549
## [25] -0.36235701 -0.28518413 -0.28215348 -0.24353279 -0.17108968 -0.10703609
## [31] -0.10492484 -0.06990653 0.00000000 0.00000000 0.00000000 -0.49553611
Por otro lado, la expresión del modelo como MA infinito tiene por coeficientes los siguientes:
##Pesos Psi (MA infinit)
(psis=ARMAtoMA(ar=mod2$model$phi,ma=mod2$model$theta,lag.max=36))
## [1] -0.5786551167 -0.1205744686 -0.1172765473 -0.0053997140 0.0626860555
## [6] 0.0530893321 -0.0469317727 0.0148084533 0.0708273651 -0.0303219421
## [11] -0.0143531970 -0.8130415863 0.4512073900 0.1074930490 0.1031669144
## [16] -0.0001781322 -0.0490266932 -0.0383539364 0.0358195472 -0.0124998519
## [21] -0.0586989305 0.0225746189 0.0127431220 0.0186587624 0.0052483013
## [26] -0.0096015939 -0.0082359316 0.0034078336 -0.0004417101 -0.0031299279
## [31] 0.0008570483 0.0007322081 0.0022676351 0.0011754800 -0.0010970761
## [36] -0.0012125942
Igual que para \(mod1\), se deberán comprobar las sumas de los coeficientes, tanto de sus cuadrados como de sus valores absolutos, para poder estudiar si este es estacionario y/o invertible.
sum(pis^2)
## [1] 3.364806
sum(abs(pis))
## [1] 8.191142
sum(psis^2)
## [1] 1.275297
sum(abs(psis))
## [1] 2.875794
Los cuatro sumatorios convergen a números finitos. Por tanto, la conclusión es que este segundo modelo, al igual que el primero, es estacionario e invertible, tanto expresado como un AR infinito, como expresado como un MA infinito.
En cuanto a las medidas de adecuación, se analizarán de nuevo los criterios AIC y BIC:
print(c("AIC del modelo 2:",mod2$aic))
## [1] "AIC del modelo 2:" "-1366.37233772363"
print(c("BIC del modelo 2:",BIC(mod2)))
## [1] "BIC del modelo 2:" "-1328.26092779386"
A modo de comparativa, se ha observado que ambos modelos son invertibles y causales, sin embargo las medidas de adecuación son ligeramente distintas: - mod1: AIC = -1370.72, BIC = -1355.47 - mod2: AIC = -1366.37, BIC = -1328.26
Se eliminan primero de todo las 12 últimas observaciones de la serie con la finalidad de estudiar y comparar la capacidad de predicción de ambos modelos.
ultim=c(2017,12);
ser2=window(ser, end=ultim)
Una vez definida la nueva serie se ajustan ambos modelos a ésta.
lnser2=log(ser2)
d1d12lnser2=diff(diff(lnser2,12))
(mod1_1 <- arima(d1d12lnser2, order=c(0,0,3), seasonal=list(order=c(0,0,1),period=12)))
##
## Call:
## arima(x = d1d12lnser2, order = c(0, 0, 3), seasonal = list(order = c(0, 0, 1),
## period = 12))
##
## Coefficients:
## ma1 ma2 ma3 sma1 intercept
## -0.5688 -0.1065 -0.0832 -0.8046 -1e-04
## s.e. 0.0566 0.0652 0.0587 0.0406 1e-04
##
## sigma^2 estimated as 0.000914: log likelihood = 665.12, aic = -1318.23
(mod1_1 <- arima(lnser2, order=c(0,1,3), seasonal=list(order=c(0,1,1),period=12)))
##
## Call:
## arima(x = lnser2, order = c(0, 1, 3), seasonal = list(order = c(0, 1, 1), period = 12))
##
## Coefficients:
## ma1 ma2 ma3 sma1
## -0.5681 -0.1054 -0.0813 -0.7974
## s.e. 0.0567 0.0651 0.0585 0.0399
##
## sigma^2 estimated as 0.0009172: log likelihood = 664.76, aic = -1319.52
cat("\nT-ratios", round(mod1_1$coef/sqrt(diag(mod1_1$var.coef)),2))
##
## T-ratios -10.02 -1.62 -1.39 -19.97
(mod2_1 <- arima(d1d12lnser2, order=c(8,0,0), seasonal=list(order=c(0,0,1),period=12)))
##
## Call:
## arima(x = d1d12lnser2, order = c(8, 0, 0), seasonal = list(order = c(0, 0, 1),
## period = 12))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.5780 -0.4456 -0.4398 -0.3876 -0.2703 -0.1639 -0.1529 -0.1028
## s.e. 0.0565 0.0646 0.0678 0.0702 0.0707 0.0681 0.0646 0.0567
## sma1 intercept
## -0.8018 -1e-04
## s.e. 0.0418 1e-04
##
## sigma^2 estimated as 0.0009006: log likelihood = 667.51, aic = -1313.01
(mod2_1 <- arima(lnser2, order=c(8,1,0), seasonal=list(order=c(0,1,1),period=12)))
##
## Call:
## arima(x = lnser2, order = c(8, 1, 0), seasonal = list(order = c(0, 1, 1), period = 12))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.5778 -0.4448 -0.4381 -0.3852 -0.2687 -0.1628 -0.1522 -0.1016
## s.e. 0.0566 0.0648 0.0679 0.0703 0.0709 0.0683 0.0647 0.0567
## sma1
## -0.7963
## s.e. 0.0411
##
## sigma^2 estimated as 0.0009029: log likelihood = 667.25, aic = -1314.5
cat("\nT-ratios", round(mod2_1$coef/sqrt(diag(mod2_1$var.coef)),2))
##
## T-ratios -10.21 -6.87 -6.45 -5.48 -3.79 -2.39 -2.35 -1.79 -19.35
Escogemos los modelos estimados a partir de las diferenciaciones hechas manualmente y eliminamos aquellos coeficientes que tengan un valor inferior a 2 en valor absoluto en el test de t-ratios. Los valores dudosos, cercanos a este valor, se comprueba que el AIC del modelo aumenta en eliminarlos y, por tanto, no deben ser fijados a 0. Una vez fijado el único parámetro que lo requiere, se comprueba que los T-ratios están todos por encima de 2, en valor absoluto, y no se han visto afectados después de ello.
(mod1_1 <- arima(lnser2, order=c(0,1,3), seasonal=list(order=c(0,1,1),period=12),fixed=c(NA,NA,0,NA)))
##
## Call:
## arima(x = lnser2, order = c(0, 1, 3), seasonal = list(order = c(0, 1, 1), period = 12),
## fixed = c(NA, NA, 0, NA))
##
## Coefficients:
## ma1 ma2 ma3 sma1
## -0.5868 -0.1588 0 -0.7951
## s.e. 0.0532 0.0498 0 0.0400
##
## sigma^2 estimated as 0.0009231: log likelihood = 663.8, aic = -1319.6
cat("\nT-ratios", round(mod1_1$coef/sqrt(diag(mod1_1$var.coef)),2))
## Warning in mod1_1$coef/sqrt(diag(mod1_1$var.coef)): longer object length is not
## a multiple of shorter object length
##
## T-ratios -11.03 -3.19 0 -14.95
rbind(coef(mod1),coef(mod1_1))
## ma1 ma2 ma3 sma1
## [1,] -0.5861043 -0.1635852 0 -0.7923814
## [2,] -0.5867856 -0.1587501 0 -0.7950725
rbind(coef(mod2),coef(mod2_1))
## ar1 ar2 ar3 ar4 ar5 ar6
## [1,] -0.5786551 -0.4554162 -0.4505765 -0.3889024 -0.2732165 -0.1709281
## [2,] -0.5777656 -0.4447952 -0.4381295 -0.3851880 -0.2687443 -0.1628127
## ar7 ar8 sma1
## [1,] -0.1675566 -0.1116351 -0.7913315
## [2,] -0.1522226 -0.1016473 -0.7963182
Se concluye que, pese a eliminar las 12 últimas observaciones, los coeficientes de ambos modelos tienen tanto el mismo signo como un valor similar, de manera que son estables y podemos confiar en las predicciones de estos, ya que son similares a los originales.
A continuación, a través de la función \(accuracy()\) de la librería \(forecast\) se evaluará la capacidad de previsión de ambos modelos, generando una predicción para los próximos 12 períodos a partir de los modelos ajustados con la serie recortada, y evaluando su precisión respecto a las observaciones originales. La salida es una tabla que contiene diferentes errores de predicción, tales como el MAPE o el MAE. Cuanto menor sean, mayor capacidad de previsión del modelo.
# Generar predicciones del modelo para los próximos 12 períodos
forecast_mod1 <- forecast(mod1_1, h = 12)
# Evaluar la precisión de las predicciones respecto a las observaciones originales
accuracy(forecast_mod1, lnser)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001416189 0.0298461 0.02289216 -0.01465075 0.2326499 0.5502584
## Test set -0.004921900 0.0211000 0.01778297 -0.04891289 0.1761067 0.4274488
## ACF1 Theil's U
## Training set 0.003397685 NA
## Test set 0.241318167 0.353361
# Generar predicciones del modelo para los próximos 12 períodos
forecast_mod2 <- forecast(mod2_1, h = 12)
# Evaluar la precisión de las predicciones respecto a las observaciones originales
accuracy(forecast_mod2, lnser)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001323747 0.02951843 0.02264060 -0.01374649 0.2301263 0.5442115
## Test set -0.007239499 0.02083684 0.01844969 -0.07190110 0.1827671 0.4434748
## ACF1 Theil's U
## Training set -0.006918027 NA
## Test set 0.229752766 0.3487475
A través de las tablas se puede observar como todas las medidas de error tienen valores similares para los dos modelos, siendo en algunas muy ligeramente superior el primer modelo, así como a la inversa en otras. De todos modos, las medidas de error tienen valores bajos para ambos modelos, así que se puede afirmar que ambos tendrán buena capacidad de previsión. Teniendo en cuenta también que son valores muy similares, se puede concluir que no será un factor determinante a la hora de escoger modelo, pues sea cual sea el elegido, la bondad de precisión no variará apenas de uno a otro.
\(\textbf{Selección de modelo}\)
El último paso en la validación de modelos es escoger con cuál de los dos ajustados se realizarán las predicciones. En primer lugar, respecto al análisis de residuos, se ha podido observar cómo ambos tienen las mismas características: varianza constante, independencia y ruido blanco y, en ninguno de ellos se cumple la hipótesis de normalidad de residuos. Respecto al estudio de la invertibilidad, causalidad y estabilidad de los modelos, también se ha podido concluir que comparten estas características. Ambos son invertibles, estables y causales. Además de todo el parecido en estas características, se ha podido ver en el último análisis como su capacidad de predicción es muy similar. Por tanto, ninguno de estos rasgos será completamente determinante a la hora de decidir entre dos modelos. En lo que sí hay diferencia entre ambos es en las medidas de adecuación. No es una diferencia abismal, pero teniendo en cuenta la importancia de estas, y que en todas las demás facetas son prácticamente iguales, será diferencial a la hora de escoger.
(AIC1=AIC(mod1))
## [1] -1370.715
(BIC1=BIC(mod1))
## [1] -1355.471
Por tanto, el modelo escogido para llevar a cabo las predicciones será el primer modelo, \(mod1\): \[ARIMA(0,1,3)(0,1,1)_{12}\]
Para el modelo ajustado sin las últimas 12 observaciones, obtenemos las predicciones puntuales, y el correspondiente intervalo de confianza al 95% para el último año. En ambas predicciones se grafican enventanadas desde el año 2015, de manera que pueda verse clara tanto la predicción como su intervalo de confianza. Graficando la serie completa no se ve de manera clara, pues los datos van desde el 1990 hasta el año que se está prediciendo.
Sabemos que el intervalo de confianza del 95% corresponderá al valor predicho \(\pm\) 1.96, que es el valor de la distribución normal que contiene el 95%, multiplicado por el error estándar. \[x_{pred} \pm 1.96 se\]
Como en el modelo se ha ajustado en una escala logarítmica, de cara a predecir será necesario invertir la escala en la que se trabaja, de manera que con la función exponencial se recuperarán los valores originales de la serie .
En primer lugar, las predicciones con la serie recortada sobre el último año observado en la serie original son las siguientes:
(p1_1<-predict(mod1_1,n.ahead=12))
## $pred
## Jan Feb Mar Apr May Jun Jul Aug
## 2018 10.18074 10.09249 10.11777 10.01816 10.04227 10.07729 10.16220 10.11308
## Sep Oct Nov Dec
## 2018 10.07201 10.05958 10.08129 10.12994
##
## $se
## Jan Feb Mar Apr May Jun
## 2018 0.03038262 0.03287431 0.03377118 0.03464485 0.03549702 0.03632921
## Jul Aug Sep Oct Nov Dec
## 2018 0.03714276 0.03793886 0.03871860 0.03948295 0.04023277 0.04096888
(pr<-exp(p1_1$pred))
## Jan Feb Mar Apr May Jun Jul Aug
## 2018 26389.87 24160.87 24779.33 22430.15 22977.48 23796.39 25905.26 24663.41
## Sep Oct Nov Dec
## 2018 23671.12 23378.79 23891.69 25082.82
ul <- exp(p1_1$pred+1.96*p1_1$se)
ll <- exp(p1_1$pred-1.96*p1_1$se)
ts.plot(ser,ll,ul,pr,lty=c(1,2,2,1),col=c(1,4,4,2),xlim=c(2015,2020), type="o"); abline(v=2012:2020,lty=3,col=4)
Se puede observar como la predicción, graficada en rojo, se encuentra
dentro del intervalo de confianza, graficado en azul, y se ajusta de una
manera precisa a las observaciones originales, graficadas en negro.
Mientras que la predicción sobre la serie original a un año vista es la siguiente:
(p1<-predict(mod1,n.ahead=12))
## $pred
## Jan Feb Mar Apr May Jun Jul Aug
## 2018
## 2019 10.17943 10.09706 10.12935 10.02755 10.04526 10.07474 10.16144 10.12179
## Sep Oct Nov Dec
## 2018 10.12319
## 2019 10.07463 10.06367 10.08450
##
## $se
## Jan Feb Mar Apr May Jun
## 2018
## 2019 0.03262028 0.03348137 0.03432087 0.03514031 0.03594108 0.03672439
## Jul Aug Sep Oct Nov Dec
## 2018 0.03014060
## 2019 0.03749134 0.03824291 0.03897999 0.03970339 0.04041385
(pr<-exp(p1$pred))
## Jan Feb Mar Apr May Jun Jul Aug
## 2018
## 2019 26355.40 24271.60 25067.94 22641.68 23046.30 23735.79 25885.58 24879.19
## Sep Oct Nov Dec
## 2018 24914.13
## 2019 23733.15 23474.43 23968.51
ul <- exp(p1$pred+1.96*p1$se)
ll <- exp(p1$pred-1.96*p1$se)
ts.plot(ser,ll,ul,pr,lty=c(1,2,2,1),col=c(1,4,4,2),xlim=c(2015,2020), type="o"); abline(v=2012:2020,lty=3,col=4)
A continuación se obtienen las medidas de error RMSE, MAE, RMSPE y MAPE
del modelo. Alguna de ellas ya ha sido mostrada anteriormente en la
evaluación de la capacidad de previsión del modelo.
obs=window(ser,start=2017)
pr<-exp(p1_1$pred)
(RMSE1=sqrt(mean((obs-pr)^2)))
## [1] 517.8877
(MAE1=mean(abs(obs-pr)))
## [1] 434.0531
(RMSPE1=sqrt(mean(((obs-pr)/obs)^2)))
## [1] 0.02119971
(MAPE1=mean(abs(obs-pr)/obs))
## [1] 0.0178454
(CI1=mean(ul-ll))
## [1] 3441.743
En las series temporales hay ciertas configuraciones en el mes que pueden afectar las predicciones, a eso se le conoce como efectos de calendario. Se considerarán dos casos principales: Semana Santa y la configuración de los días de la semana.
El primero porque dependiendo del año, la Semana Santa puede caer en el mes de Marzo, en el mes de Abril, o en ambos, y eso podría no verse reflejado en las predicciones y que estas fallaran. Con tal de lidiar con este problema, se intentará cambiar la serie con tal de que la Semana Santa sea considerada la mitad de días para cada mes.
El segundo porque según la proporcion de días laborables y fines de semana a lo largo de los meses de un año puede afectar también a las predicciones hechas sobre la serie. En este caso, se lidiará con el problema estableciendo la proporción de Trading Days/Weekends a 5/2 en todos los meses.
source("CalendarEffects.r")
inici=c(1990,1,length(ser))
(vEa=Weaster(inici))
## Jan Feb Mar Apr May Jun
## 1990 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 1991 0.0000000 0.0000000 0.5000000 -0.5000000 0.0000000 0.0000000
## 1992 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 1993 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 1994 0.0000000 0.0000000 0.1666667 -0.1666667 0.0000000 0.0000000
## 1995 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 1996 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 1997 0.0000000 0.0000000 0.5000000 -0.5000000 0.0000000 0.0000000
## 1998 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 1999 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2000 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2001 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2002 0.0000000 0.0000000 0.5000000 -0.5000000 0.0000000 0.0000000
## 2003 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2004 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2005 0.0000000 0.0000000 0.5000000 -0.5000000 0.0000000 0.0000000
## 2006 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2007 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2008 0.0000000 0.0000000 0.5000000 -0.5000000 0.0000000 0.0000000
## 2009 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2010 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2011 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2012 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2013 0.0000000 0.0000000 0.5000000 -0.5000000 0.0000000 0.0000000
## 2014 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2015 0.0000000 0.0000000 -0.1666667 0.1666667 0.0000000 0.0000000
## 2016 0.0000000 0.0000000 0.5000000 -0.5000000 0.0000000 0.0000000
## 2017 0.0000000 0.0000000 -0.5000000 0.5000000 0.0000000 0.0000000
## 2018 0.0000000 0.0000000 0.5000000 -0.5000000 0.0000000 0.0000000
## Jul Aug Sep Oct Nov Dec
## 1990 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1991 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1992 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1993 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1994 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1995 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1996 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1997 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1998 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 1999 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2001 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2002 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2003 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2004 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2005 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2006 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2007 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2008 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2009 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2010 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2011 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2012 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2013 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2014 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2015 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2016 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2017 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 2018 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
En este listado se puede ver la repartición de los días de la Semana Santa entre Marzo y Abril para cada año. Se puede ver como, por ejemplo, en el año 1993 Abril posee el doble de días que Marzo para estas fechas (Marzo -0.5, Abril +0.5). Esto quiere decir que la Semana Santa cayó en Abril entera, por lo que en ese mes hay el doble de días respecto a la repartición deseada. En cambio, se puede ver como en 1994 la repartición de días fue 5 días de Semana Santa en Marzo y 3 en Abril, es decir un 17% más en el primer mes de los dos, respecto a la repartición deseada. Los años en que la proporción es $$0, es porque están igualmente repartidos estos días entre los dos meses.
(vTD=Wtrad(inici))
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1990 3.0 0.0 -0.5 -1.5 3.0 -1.5 -0.5 3.0 -5.0 3.0 2.0 -4.0
## 1991 3.0 0.0 -4.0 2.0 3.0 -5.0 3.0 -0.5 -1.5 3.0 -1.5 -0.5
## 1992 3.0 -2.5 -0.5 2.0 -4.0 2.0 3.0 -4.0 2.0 -0.5 -1.5 3.0
## 1993 -4.0 0.0 3.0 2.0 -4.0 2.0 -0.5 -0.5 2.0 -4.0 2.0 3.0
## 1994 -4.0 0.0 3.0 -1.5 -0.5 2.0 -4.0 3.0 2.0 -4.0 2.0 -0.5
## 1995 -0.5 0.0 3.0 -5.0 3.0 2.0 -4.0 3.0 -1.5 -0.5 2.0 -4.0
## 1996 3.0 1.0 -4.0 2.0 3.0 -5.0 3.0 -0.5 -1.5 3.0 -1.5 -0.5
## 1997 3.0 0.0 -4.0 2.0 -0.5 -1.5 3.0 -4.0 2.0 3.0 -5.0 3.0
## 1998 -0.5 0.0 -0.5 2.0 -4.0 2.0 3.0 -4.0 2.0 -0.5 -1.5 3.0
## 1999 -4.0 0.0 3.0 2.0 -4.0 2.0 -0.5 -0.5 2.0 -4.0 2.0 3.0
## 2000 -4.0 1.0 3.0 -5.0 3.0 2.0 -4.0 3.0 -1.5 -0.5 2.0 -4.0
## 2001 3.0 0.0 -0.5 -1.5 3.0 -1.5 -0.5 3.0 -5.0 3.0 2.0 -4.0
## 2002 3.0 0.0 -4.0 2.0 3.0 -5.0 3.0 -0.5 -1.5 3.0 -1.5 -0.5
## 2003 3.0 0.0 -4.0 2.0 -0.5 -1.5 3.0 -4.0 2.0 3.0 -5.0 3.0
## 2004 -0.5 -2.5 3.0 2.0 -4.0 2.0 -0.5 -0.5 2.0 -4.0 2.0 3.0
## 2005 -4.0 0.0 3.0 -1.5 -0.5 2.0 -4.0 3.0 2.0 -4.0 2.0 -0.5
## 2006 -0.5 0.0 3.0 -5.0 3.0 2.0 -4.0 3.0 -1.5 -0.5 2.0 -4.0
## 2007 3.0 0.0 -0.5 -1.5 3.0 -1.5 -0.5 3.0 -5.0 3.0 2.0 -4.0
## 2008 3.0 1.0 -4.0 2.0 -0.5 -1.5 3.0 -4.0 2.0 3.0 -5.0 3.0
## 2009 -0.5 0.0 -0.5 2.0 -4.0 2.0 3.0 -4.0 2.0 -0.5 -1.5 3.0
## 2010 -4.0 0.0 3.0 2.0 -4.0 2.0 -0.5 -0.5 2.0 -4.0 2.0 3.0
## 2011 -4.0 0.0 3.0 -1.5 -0.5 2.0 -4.0 3.0 2.0 -4.0 2.0 -0.5
## 2012 -0.5 1.0 -0.5 -1.5 3.0 -1.5 -0.5 3.0 -5.0 3.0 2.0 -4.0
## 2013 3.0 0.0 -4.0 2.0 3.0 -5.0 3.0 -0.5 -1.5 3.0 -1.5 -0.5
## 2014 3.0 0.0 -4.0 2.0 -0.5 -1.5 3.0 -4.0 2.0 3.0 -5.0 3.0
## 2015 -0.5 0.0 -0.5 2.0 -4.0 2.0 3.0 -4.0 2.0 -0.5 -1.5 3.0
## 2016 -4.0 1.0 3.0 -1.5 -0.5 2.0 -4.0 3.0 2.0 -4.0 2.0 -0.5
## 2017 -0.5 0.0 3.0 -5.0 3.0 2.0 -4.0 3.0 -1.5 -0.5 2.0 -4.0
## 2018 3.0 0.0 -0.5 -1.5 3.0 -1.5 -0.5 3.0 -5.0 3.0 2.0
En el listado para los Trading Days se puede observar cuantos días laborables hay de más, o de menos, respecto a la proporción de 5/2 buscada. Se puede ver como en 1990, por ejemplo, hubo 3 días laborables más de los que tocarían respecto a esta partición, así como en 1993 hubo 4 días más de fin de semana respecto a ella.
A continuación, se ajustará un modelo para cada efecto de calendario, además de otro con ambos efectos juntos. De esta manera, se podrá estudiar su efecto respecto al modelo principal ajustado anteriormente:
(mod1=arima(lnser,order=c(0,1,2),
seasonal=list(order=c(0,1,1),period=12)))
##
## Call:
## arima(x = lnser, order = c(0, 1, 2), seasonal = list(order = c(0, 1, 1), period = 12))
##
## Coefficients:
## ma1 ma2 sma1
## -0.5861 -0.1636 -0.7924
## s.e. 0.0522 0.0489 0.0383
##
## sigma^2 estimated as 0.0009085: log likelihood = 689.36, aic = -1370.72
(mod1Ea=arima(lnser,order=c(0,1,2),
seasonal=list(order=c(0,1,1),period=12),
xreg=data.frame(vEa)))
##
## Call:
## arima(x = lnser, order = c(0, 1, 2), seasonal = list(order = c(0, 1, 1), period = 12),
## xreg = data.frame(vEa))
##
## Coefficients:
## ma1 ma2 sma1 vEa
## -0.5659 -0.1794 -0.7863 -0.0222
## s.e. 0.0515 0.0482 0.0386 0.0063
##
## sigma^2 estimated as 0.0008777: log likelihood = 695.27, aic = -1380.54
(mod1TD=arima(lnser,order=c(0,1,2),
seasonal=list(order=c(0,1,1),period=12),
xreg=data.frame(vTD)))
##
## Call:
## arima(x = lnser, order = c(0, 1, 2), seasonal = list(order = c(0, 1, 1), period = 12),
## xreg = data.frame(vTD))
##
## Coefficients:
## ma1 ma2 sma1 vTD
## -0.5718 -0.1727 -0.7877 0.0015
## s.e. 0.0518 0.0485 0.0390 0.0004
##
## sigma^2 estimated as 0.0008774: log likelihood = 695.3, aic = -1380.6
(mod1EC=arima(lnser,order=c(0,1,2),
seasonal=list(order=c(0,1,1),period=12),
xreg=data.frame(vEa,vTD)))
##
## Call:
## arima(x = lnser, order = c(0, 1, 2), seasonal = list(order = c(0, 1, 1), period = 12),
## xreg = data.frame(vEa, vTD))
##
## Coefficients:
## ma1 ma2 sma1 vEa vTD
## -0.5531 -0.1877 -0.7830 -0.0201 0.0014
## s.e. 0.0513 0.0480 0.0391 0.0062 0.0004
##
## sigma^2 estimated as 0.0008518: log likelihood = 700.36, aic = -1388.72
Después de ajustar los tres modelos con efectos de calendario, se puede observar como a cada cual disminuye tanto el AIC como la varianza estimada. En este caso, el efecto causado por los Trading Days es prácticamente idéntico al causado por la Semana Santa, pero el efecto provocado por los dos juntos es aun mayor que los dos anteriores por separado. En este último modelo se puede ver como los ratios de las variables que contienen estos efectos (-3.23, 3.20) són mayores que 2, en valor absoluto, y, por tanto, significativos.
cat("\nT-ratios", round(mod1EC$coef/sqrt(diag(mod1EC$var.coef)),2))
##
## T-ratios -10.79 -3.91 -20.03 -3.23 3.2
modelo ajustado para la Semana Santa: var. est: 0.0008777 ; AIC: -1380.54 modelo ajustado para los Trading Days: var. est: 0.0008774 ; AIC: -1380.6 modelo ajustado para los dos: var. est: 0.0008518 ; AIC: -1388.72
Una vez estudiados los efectos de calendario sobre el modelo, es conveniente estudiar también el efecto de alguna circunstancia concreta a lo largo de los años que comprende la serie que pudiera haber causado cambios en esta.
Entre 1990 y 2018 destacan dos factores que pudieran tener efecto sobre la misma: - La liberalización del mercado eléctrico de 1998, en búsqueda de una mayor competencia. - La introducción del Plan de Acción de Energías Renovables en 2005, el cual se introdujo con objetivo de aumentar la proporción de energías renovables en el consumo eléctrico del país.
Se ajustará un modelo añadiendo el efecto de estos, de tal manera que pueda verse si el efecto causado en la serie es suficientemente significativo como para ser considerado en este.
vLib=ts(rep(0,length(ser)),start=1990,freq=12)
window(vLib,start=c(1998,1))<-1
# Inicio en enero de 1998 con la libre elegibilidad para los grandes consumidores de energía
vPlan=ts(rep(0,length(ser)),start=1990,freq=12)
window(vPlan,start=c(2005,8))<-1
# Plan aprobado en agosto(8) de 2005
(mod1ECIA=arima(lnser,order=c(0,1,2),
seasonal=list(order=c(0,1,1),period=12),
xreg=data.frame(vEa,vTD,vLib,vPlan)))
##
## Call:
## arima(x = lnser, order = c(0, 1, 2), seasonal = list(order = c(0, 1, 1), period = 12),
## xreg = data.frame(vEa, vTD, vLib, vPlan))
##
## Coefficients:
## ma1 ma2 sma1 vEa vTD vLib vPlan
## -0.5462 -0.1777 -0.7849 -0.0201 0.0014 -0.0152 -0.0243
## s.e. 0.0523 0.0490 0.0387 0.0062 0.0004 0.0215 0.0209
##
## sigma^2 estimated as 0.0008472: log likelihood = 701.26, aic = -1386.52
cat("\nT-ratios", round(mod1ECIA$coef/sqrt(diag(mod1ECIA$var.coef)),2))
##
## T-ratios -10.44 -3.63 -20.29 -3.24 3.25 -0.71 -1.16
Después del ajuste realizado, se puede ver como, por un lado, el AIC de este último modelo es mayor que el del modelo que contiene solo efectos de calendario. Mientras que, por otro lado, al mostrar los ratios entre sus coeficientes y su error estándar, se concluye que no son parámetros influyentes, pues estos son menores que 2, en valor absoluto. Pese a que después de un análisis exhaustivo de circunstancias que pudieran afectar al consumo eléctrico en España se hubieran encontrado los mencionados como factores principales, se concluye que estos no son significativamente distintos de 0 y, por tanto, no deben ser incluidos en el modelo. Su efecto asociado a la serie temporal no es diferencial a la hora de realizar buenas predicciones.
Una vez queda seleccionado qué modelo se ajusta mejor a todos estos fenómenos, se realizará una comparativa de los valores de la serie original con los valores de la serie que contiene estos efectos.
exp(coef(mod1EC)["vEa"])
## vEa
## 0.9800729
exp(coef(mod1EC)["vTD"])
## vTD
## 1.00137
lnserEC=lnser-
coef(mod1EC)["vEa"]*vEa-
coef(mod1EC)["vTD"]*vTD
plot(lnser, xlab="Tiempo", ylab="GwH(log)", main="Comaparativa del logaritmo de la serie con y sin efectos de calendario aplicados")
lines(lnserEC, col="red")
abline(v=1985:2019,lty=3,col=4)
Una vez vista la comparativa, y visto que los efectos de calendario no
cambian la serie, deberemos repetir todo el proceso de identificación,
estimación y validación de un modelo que ajuste la nueva serie
En primer lugar se han de realizar las mismas diferenciaciones que con la serie, pues esta se comportará igual ante tales transformaciones. Es decir, una regular y una estacional.
d1d12lnserEC=diff(diff(lnserEC,12))
par(mfrow=c(1,2))
acf(d1d12lnserEC,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
pacf(d1d12lnserEC,ylim=c(-1,1),lag.max=72,col=c(rep(1,11),2),lwd=2)
par(mfrow=c(1,1))
De la representación del ACF y PACF de la serie se extraerán los parámetros de los modelos candidatos y, al igual que con la serie original, se ajustarán dos para decidir cuál es el que mejor se ajusta a los datos.
Parte estacional: Del ACF se pueden extraer dos retardos significativos. Por tanto, se planteará un MA(1) estacional. Por otro lado, para el PACF se observan también dos, el último que aparece en la representación se considera que es un valor fruto del azar. Pot tanto, se planteará un AR(2) estacional.
Parte regular: En el ACF se puede observar que el último retardo significativo es el tercero. Los posteriores, como antes, se considerarán fruto del azar. Se considerará un MA(3) regular. Por otro lado, para el PACF se pueden observar 8 retardos significativos, por lo que el modelo será un AR(8).
Tal como se ha comentado, se observan los mismos retardos significativos tanto para la parte regular como para la parte estacional que en los ACF y PACF de la serie estacionaria original. Por ende, los modelos a estimar serán los mismos. Puesto que el modelo seleccionado anteriormente ha sido el que contiene un MA(2) regular y un MA(1) estacional, cabe esperar que este sea el de mejor ajuste. De todas maneras, se probará a estimar un \(ARIMA(8,1,0)(0,1,1)_{12}\) con los efectos de calendario, para comprobar que el AIC del modelo es mayor y se ha de proseguir el análisis con el modelo ya ajustado
(modECcand=arima(lnser,order=c(8,1,0),
seasonal=list(order=c(0,1,1),period=12),
xreg=data.frame(vEa,vTD)))
##
## Call:
## arima(x = lnser, order = c(8, 1, 0), seasonal = list(order = c(0, 1, 1), period = 12),
## xreg = data.frame(vEa, vTD))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8
## -0.5371 -0.4357 -0.4470 -0.3515 -0.2734 -0.1802 -0.1378 -0.1107
## s.e. 0.0553 0.0625 0.0657 0.0682 0.0690 0.0659 0.0628 0.0558
## sma1 vEa vTD
## -0.7819 -0.0206 0.0013
## s.e. 0.0405 0.0062 0.0004
##
## sigma^2 estimated as 0.0008327: log likelihood = 704.07, aic = -1384.14
cat("\nT-ratios", round(modECcand$coef/sqrt(diag(modECcand$var.coef)),2))
##
## T-ratios -9.71 -6.97 -6.81 -5.16 -3.96 -2.73 -2.2 -1.98 -19.3 -3.32 3.11
modEC = mod1EC
Efectivamente, el AIC de este modelo es mayor que del ya ajustado previamente, que tenía AIC = -1388.72.
Pese a ser el mismo modelo, se ha de analizar los residuos y comprobar que se cumplen las hipótesis de normalidad, varianza constante e independencia.
Antes de comprobar estas hipótesis, se comprobara a través del ACF y PACF que estos caigan dentro de las bandas de confianza, como se espera que lo hagan.
resi=resid(modEC)
par(mfrow=c(1,2))
acf(resi,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
pacf(resi,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
par(mfrow=c(1,1))
plot(resi)
abline(h=0)
abline(h=c(-3*sd(resi),3*sd(resi)),lty=3,col=c(1))
Tal como se esperaba todas la observaciones, o prácticamente todas, caen
dentro de las bandas de confianza. Las que sobresalen por la mínima, se
consideran valores fruto del azar.
En la primera exploración se observa que la variabilidad e similar para la serie, sin contar el pico pronunciado previo al 1995. Por tanto, se espera que la varianza sea constante. Para reforzar esta hipótesis se realizará un scatterplot.
scatter.smooth(sqrt(abs(resi)),lpars=list(col=2))
Los residuos se distribuyen homogéneamente alrededor de la línea roja,
que es prácticamente horizontal. Se observan algunas observaciones
atípicas, que serán tratadas en el siguiente apartado. La varianza de
los residuos es constante.
A continuación, se comprobará la normalidad de los residuos. En primer lugar, con el Normal Q-Q Plot.
qqnorm(resi)
qqline(resi,col=2,lwd=2)
No hay volatilidad, las colas en los extremos son de observaciones
atípicas. Los valores se ajustan en una recta, por lo que de momento se
cumple la hipótesis. El test que es determinante para concluir esto, es
el de Shapiro-Wilks.
hist(resi,breaks=20,freq=FALSE)
curve(dnorm(x,mean=mean(resi),sd=sd(resi)),col=2,add=T)
shapiro.test(resi)
##
## Shapiro-Wilk normality test
##
## data: resi
## W = 0.99022, p-value = 0.02065
El histograma es solo ilustrativo, el test da un p-valor inferior a 0.05, por lo que se rechaza la hipótesis de normalidad para los residuos del modelo.
A través del test de Llung-Box se puede comprobar la independencia de estos, tal como se ha visto anteriormente.
tsdiag(modEC,gof.lag=72)
Las observaciones son independientes y el ruido es blanco, todas se
sitúan por encima de la línea discontínua azul. De este modo se concluye
el análisis de los residuos, pudiendo afirmar que el modelo estimado
explica suficientemente bien los datos.
El modelo será causal si y solo si, los módulos de las raíces del polinomio característico respecto a las phis son mayores que 1 (>1)
Mod(polyroot(c(1,-modEC$model$phi)))
## numeric(0)
Como era de esperar, pues ya ha pasado con el modelo ajustado sobre la serie original transformada, no tiene polinomio característico respecto a la parte autoregresiva del mismo, al tratarse de un MA regular. Por tanto, el modelo es automáticamente causal. Por otra banda, el modelo será invertible si el módulo de las raíces del polinomio característico respecto a las thetas son mayores que 1 (>1)
Mod(polyroot(c(1,modEC$model$theta)))
## [1] 1.020597 1.020597 1.020597 1.020597 1.020597 1.020597 1.020597 1.020597
## [9] 1.020597 1.020597 1.020597 1.020597 1.264989 4.211103
Efectivamente, es el caso y, por tanto, el modelo es también invertible. En cuanto a las medidas de adecuación del mismo, son las siguientes.
(AIC2=AIC(modEC))
## [1] -1388.719
(BIC2=BIC(modEC))
## [1] -1365.852
Tanto el AIC como el BIC de este modelo son menores que los del modelo sin efectos de calendario estimado anteriormente.
Como se está comprobando, el modelo tiene las mismas propiedades tanto para la serie original como para la que incluye efectos de calendario. Por consecuencia, también compartirá la estabilidad y la capacidad de previsión al enventanar la serie sin la última observación. De modo que se pasará directamente a las predicciones, tal como se ha hecho anteriormente.
En primer lugar se retiran de la serie las últimas 12 observaciones, y se hace también un enventanado de las componentes que contienen los efectos de calendario, con tal de modelar también incluyéndolas.
ultim=c(2017,12)
ser2=window(ser, end=ultim)
lnser2=log(ser2)
vEa2=window(vEa,end=ultim)
vTD2=window(vTD,end=ultim)
modEC
##
## Call:
## arima(x = lnser, order = c(0, 1, 2), seasonal = list(order = c(0, 1, 1), period = 12),
## xreg = data.frame(vEa, vTD))
##
## Coefficients:
## ma1 ma2 sma1 vEa vTD
## -0.5531 -0.1877 -0.7830 -0.0201 0.0014
## s.e. 0.0513 0.0480 0.0391 0.0062 0.0004
##
## sigma^2 estimated as 0.0008518: log likelihood = 700.36, aic = -1388.72
(modEC2=arima(lnser2,order=c(0,1,3),
seasonal=list(order=c(0,1,1),period=12),
xreg=data.frame(vEa2,vTD2)))
##
## Call:
## arima(x = lnser2, order = c(0, 1, 3), seasonal = list(order = c(0, 1, 1), period = 12),
## xreg = data.frame(vEa2, vTD2))
##
## Coefficients:
## ma1 ma2 ma3 sma1 vEa2 vTD2
## -0.5242 -0.1236 -0.1002 -0.7862 -0.0222 0.0013
## s.e. 0.0569 0.0614 0.0564 0.0413 0.0064 0.0004
##
## sigma^2 estimated as 0.000853: log likelihood = 676.77, aic = -1339.54
vEa2=window(vEa,start=ultim+c(0,1))
vTD2=window(vTD,start=ultim+c(0,1))
A continuación se realizarán las predicciones, con sus respectivos intervalos de confianza, tanto con el modelo original (entiéndase por original el modelo sin recortar las últimas observaciones) como con el último estimado.
pre=predict(modEC2,n.ahead = 12,newxreg=data.frame(vEa2,vTD2))
## Warning in z[[1L]] + xm: longer object length is not a multiple of shorter
## object length
ll=exp(pre$pred-1.96*pre$se)
pr=exp(pre$pred)
ul=exp(pre$pred+1.96*pre$se)
ts.plot(ser,ll,ul,pr,
lty=c(1,2,2,1),
col=c(1,4,4,2),
xlim=c(2016,2020), type="o")
abline(v=2016:2020,lty=3,col=4)
La predicción de la serie modelada con \(modEC2\), representada en rojo, se
encuentra dentro de su intervalo de confianza, representado en azul, y
se ajusta de una manera muy precisa al gráfico de la serie original,
representado en negro.
pre=predict(modEC,n.ahead = 12,newxreg=data.frame(vEa2,vTD2))
## Warning in z[[1L]] + xm: longer object length is not a multiple of shorter
## object length
ll=exp(pre$pred-1.96*pre$se)
pr=exp(pre$pred)
ul=exp(pre$pred+1.96*pre$se)
ts.plot(ser,ll,ul,pr,
lty=c(1,2,2,1),
col=c(1,4,4,2),
xlim=c(2016,2020), type="o")
abline(v=2016:2020,lty=3,col=4)
En este caso la representación es de la predicción de la serie a un año
vista del final de la misma, modelada con los efectos de calendario
incluidos, es decir mediante \(modEC\).
No se puede comprobar la precisión de esta previsión, pues no tenemos
datos para ello, pero teniendo en cuenta el buen ajuste del gráfico
anterior, se puede confiar en que la predicción será precisa.
Una vez hechos los pronósticos, es necesario calcular las medidas de error del modelo, las cuales evalúan su capacidad de previsión, es decir, cuán buenas pueden ser las predicciones hechas.
pre=predict(modEC2,n.ahead = 12,newxreg=data.frame(vEa2,vTD2))
## Warning in z[[1L]] + xm: longer object length is not a multiple of shorter
## object length
pr <- exp(pre$pred)
obs=window(ser,start=2018)
(RMSE2=sqrt(mean((obs-pr)^2)))
## [1] 585.6177
(MAE2=mean(abs(obs-pr)))
## [1] 461.5682
(RMSPE2=sqrt(mean(((obs-pr)/obs)^2)))
## [1] 0.02357383
(MAPE2=mean(abs(obs-pr)/obs))
## [1] 0.01882707
(CI2=mean(ul-ll))
## [1] 3384.447
source("atipics2.r")
En primer lugar, se aplicará la detección automática de los datos atípicos.
mod.atip=outdetec(modEC,dif=c(1,12),crit=2.8,LS=T)
atipics=mod.atip$atip[order(mod.atip$atip[,1]),]
meses=c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dic")
data.frame(atipics,
Fecha=paste(meses[(atipics[,1]-1)%%12+1],start(lnser)[1]+((atipics[,1]-1)%/%12)),
PercVar=exp(atipics[,3])*100)
## Obs type_detected W_coeff ABS_L_Ratio Fecha PercVar
## 3 14 AO 0.07804566 3.436857 Feb 1991 108.11720
## 14 18 AO -0.05667709 2.889440 Jun 1991 94.48991
## 10 25 TC 0.06288620 3.041845 Ene 1992 106.49056
## 1 49 AO -0.11658634 4.947707 Ene 1994 88.99533
## 13 63 TC 0.05825792 2.925502 Mar 1995 105.99884
## 8 86 TC -0.06399438 3.015613 Feb 1997 93.80103
## 12 97 TC -0.05948991 2.951189 Ene 1998 94.22450
## 17 164 AO 0.05443965 2.874542 Ago 2003 105.59487
## 7 169 AO -0.06760858 3.154946 Ene 2004 93.46262
## 6 199 AO 0.06747490 3.104201 Jul 2006 106.98034
## 18 215 LS 0.04308703 2.861133 Nov 2007 104.40288
## 11 224 AO -0.06122103 3.008113 Ago 2008 94.06153
## 5 230 LS -0.05687798 3.207436 Feb 2009 94.47093
## 2 266 AO 0.08766443 3.785555 Feb 2012 109.16217
## 9 307 AO 0.06385094 3.019785 Jul 2015 106.59335
## 4 312 AO -0.07564278 3.314606 Dic 2015 92.71473
## 15 326 TC -0.05764212 2.878505 Feb 2017 94.39877
## 16 330 AO 0.05652395 2.851037 Jun 2017 105.81520
A través de la detección automática, se pueden detectar el tipo de outlier del que se trata (AO,TC o LS), en qué fecha se encuentran y otras medidas como:
W_coeff: Se basa en la idea de que los valores atípicos tienden a tener un impacto mayor en la varianza de una serie temporal que los valores normales. El W_coeff se calcula como la relación entre la varianza de la serie temporal que se obtiene al retirar un punto de datos específico, y la varianza de la serie temporal original. Es decir, se calcula la varianza de la serie temporal sin el punto de datos, y se divide por la varianza de la serie temporal completa. Un valor de W_coeff mayor a 1 indica que el punto de datos en cuestión tiene un impacto mayor en la varianza de la serie temporal de lo que se esperaría de un valor normal. Por lo tanto, puede ser considerado como un valor atípico potencial en la serie temporal.
ABS_L_Ratio: Se define como la relación entre el valor absoluto de la diferencia entre un punto de datos y la mediana de la serie temporal, y la distancia intercuartil (IQR) de la serie temporal. Es decir, el Abs_L ratio mide qué tan lejos se encuentra un valor de la mediana de la serie temporal en términos de la variabilidad de la serie. Un valor de Abs_L ratio mayor a 3 se considera un valor atípico potencial en una serie temporal.
PercVar: Indica la proporción de la varianza total de una serie temporal que es explicada por un modelo ajustado. Por ejemplo, si el PercVar es del 80%, significa que el 80% de la variabilidad de la serie temporal puede ser explicada por el modelo ajustado, mientras que el 20% restante representa la variabilidad que no puede ser explicada por el modelo.
Una vez detectados las observaciones atípicas, el tratamiento de estas se basa en linealizar la serie, con los efectos de calendario incluidos en ella.
lnser.lin=lineal(lnser,mod.atip$atip)
ser.lin=exp(lnser.lin)
plot(ser)
lines(exp(lnser.lin),col=2)
plot(lnser-lnser.lin)
abline(h=0, col='red',lwd = 3,lty = 2)
lnserEC.lin=lnser.lin-
coef(modEC)["vEa"]*vEa-
coef(modEC)["vTD"]*vTD
Ajuste del modelo respecto a la serie linealizada con los efectos de calendario incluidos.
d1d12lnserEC.lin=diff(diff(lnserEC.lin,12))
plot(d1d12lnserEC.lin)
abline(h=0)
par(mfrow=c(1,2))
acf(d1d12lnserEC.lin,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
pacf(d1d12lnserEC.lin,ylim=c(-1,1),lag.max=72,col=c(rep(1,11),2),lwd=2)
par(mfrow=c(1,1))
Estacional: MA(1), AR(2) Regular: MA(4), AR(5) De todas las combinaciones posibles, la que mejor se adecua a los datos, es decir, la que mejor AIC da, es la compuesta por un MA(1) estacional y un MA(4) regular
(modEClin=arima(lnser.lin,order=c(0,1,4),
seasonal=list(order=c(0,1,1),period=12),
xreg=data.frame(vEa,vTD)))
##
## Call:
## arima(x = lnser.lin, order = c(0, 1, 4), seasonal = list(order = c(0, 1, 1),
## period = 12), xreg = data.frame(vEa, vTD))
##
## Coefficients:
## ma1 ma2 ma3 ma4 sma1 vEa vTD
## -0.5720 -0.2785 -0.0313 0.1868 -0.7136 -0.0197 9e-04
## s.e. 0.0546 0.0652 0.0652 0.0555 0.0412 0.0045 3e-04
##
## sigma^2 estimated as 0.000487: log likelihood = 795.13, aic = -1574.26
cat("\nT-ratios", round(modEClin$coef/sqrt(diag(modEClin$var.coef)),2))
##
## T-ratios -10.48 -4.27 -0.48 3.37 -17.34 -4.37 3.13
(modEClin=arima(lnser.lin,order=c(0,1,4),
seasonal=list(order=c(0,1,1),period=12), fixed=c(NA,NA,0,NA,NA,NA,NA),
xreg=data.frame(vEa,vTD)))
##
## Call:
## arima(x = lnser.lin, order = c(0, 1, 4), seasonal = list(order = c(0, 1, 1),
## period = 12), xreg = data.frame(vEa, vTD), fixed = c(NA, NA, 0, NA, NA,
## NA, NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 sma1 vEa vTD
## -0.5760 -0.2914 0 0.1735 -0.714 -0.0197 1e-03
## s.e. 0.0531 0.0584 0 0.0478 0.041 0.0045 3e-04
##
## sigma^2 estimated as 0.0004873: log likelihood = 795.02, aic = -1576.03
\(\textbf{Análisis de los residuos de la serie linealizada.}\)
resi=resid(modEClin)
plot(resi)
abline(h=0)
abline(h=c(-3*sd(resi),3*sd(resi)),lty=3,col=4)
par(mfrow=c(1,2))
acf(resi,ylim=c(-1,1),lag.max=72,col=c(2,rep(1,11)),lwd=2)
pacf(resi,ylim=c(-1,1),lag.max=72,col=c(rep(1,11),2),lwd=2)
par(mfrow=c(1,1))
scatter.smooth(sqrt(abs(resi)), lpars=list(col=2))
qqnorm(resi)
qqline(resi,col=2,lwd=2)
hist(resi,breaks=15, freq=FALSE)
curve(dnorm(x, mean=mean(resi), sd=sd(resi)), col=2, lwd=2, add=T)
shapiro.test(resi)
##
## Shapiro-Wilk normality test
##
## data: resi
## W = 0.99761, p-value = 0.9009
tsdiag(modEClin,gof.lag=72)
Hipótesis de varianza constante, normalidad e independencia residual
comprobadas. El modelo queda validado y se puede utilizar para hacer
predicciones.
\(\textbf{Medidas de error}\)
(AIC3=AIC(modEClin)+2*nrow(mod.atip$atip))
## [1] -1540.034
(BIC3=BIC(modEClin)+log(length(ser)**nrow(mod.atip$atip)))
## [1] -1444.068
Más pequeños que cualquier modelo ajustado anteriormente.
\(\textbf{Causalidad e invertibilidad}\)
Mod(polyroot(c(1,-modEClin$model$phi)))
## numeric(0)
Todas las raíces por encima de 1 \(\rightarrow\) el modelo es causal.
Mod(polyroot(c(1,modEClin$model$theta)))
## [1] 1.028472 1.028472 1.028472 1.028472 1.028472 1.028472 1.028472 1.028472
## [9] 1.028472 1.028472 1.028472 1.028472 1.768972 1.768972 1.357193 1.357193
Todas las raíces por encima de 1 \(\rightarrow\) el modelo es invertible
\(\textbf{Previsiones}\)
Recortamos las últimas 12 observaciones y predecimos, tal como hemos hecho con las anteriores.
ultim=c(2017,12)
ser.lin2=window(ser.lin, end=ultim)
lnser.lin2=log(ser.lin2)
vEa2=window(vEa,end=ultim)
vTD2=window(vTD,end=ultim)
modEClin
##
## Call:
## arima(x = lnser.lin, order = c(0, 1, 4), seasonal = list(order = c(0, 1, 1),
## period = 12), xreg = data.frame(vEa, vTD), fixed = c(NA, NA, 0, NA, NA,
## NA, NA))
##
## Coefficients:
## ma1 ma2 ma3 ma4 sma1 vEa vTD
## -0.5760 -0.2914 0 0.1735 -0.714 -0.0197 1e-03
## s.e. 0.0531 0.0584 0 0.0478 0.041 0.0045 3e-04
##
## sigma^2 estimated as 0.0004873: log likelihood = 795.02, aic = -1576.03
(modEClin2=arima(lnser.lin2,order=c(6,1,0),
seasonal=list(order=c(0,1,2),period=12),
xreg=data.frame(vEa2,vTD2)))
##
## Call:
## arima(x = lnser.lin2, order = c(6, 1, 0), seasonal = list(order = c(0, 1, 2),
## period = 12), xreg = data.frame(vEa2, vTD2))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 sma1 sma2
## -0.5526 -0.5647 -0.4933 -0.3158 -0.1854 -0.0585 -0.7527 0.0451
## s.e. 0.0561 0.0635 0.0693 0.0685 0.0632 0.0571 0.0603 0.0616
## vEa2 vTD2
## -0.0212 1e-03
## s.e. 0.0048 3e-04
##
## sigma^2 estimated as 0.0004927: log likelihood = 766.81, aic = -1511.62
vEa2=window(vEa,start=ultim+c(0,1))
vTD2=window(vTD,start=ultim+c(0,1))
pre=predict(modEClin2,n.ahead = 12,newxreg=data.frame(vEa2,vTD2))
## Warning in z[[1L]] + xm: longer object length is not a multiple of shorter
## object length
wLS=sum(mod.atip$atip[mod.atip$atip[,2]=="LS",3])
ll=exp(pre$pred+wLS-1.96*pre$se)
pr=exp(pre$pred+wLS)
ul=exp(pre$pred+wLS+1.96*pre$se)
ts.plot(ser,ll,ul,pr,
lty=c(1,2,2,1),
col=c(1,4,4,2),
xlim=c(2015,2020), type="o")
abline(v=2015:2020,lty=3,col=4)
Ahora se predice con la serie original un año vista desde el final de la
serie temporal.
pre=predict(modEClin,n.ahead = 12,newxreg=data.frame(vEa2,vTD2))
## Warning in z[[1L]] + xm: longer object length is not a multiple of shorter
## object length
wLS=sum(mod.atip$atip[mod.atip$atip[,2]=="LS",3])
ll=exp(pre$pred+wLS-1.96*pre$se)
pr=exp(pre$pred+wLS)
ul=exp(pre$pred+wLS+1.96*pre$se)
ts.plot(ser,ll,ul,pr,
lty=c(1,2,2,1),
col=c(1,4,4,2),
xlim=c(2015,2020), type="o")
abline(v=2015:2020,lty=3,col=4)
Medidas de error
pre=predict(modEClin2,n.ahead = 12,newxreg=data.frame(vEa2,vTD2))
## Warning in z[[1L]] + xm: longer object length is not a multiple of shorter
## object length
pr = exp(pre$pred)
obs=window(ser,start=2018)
(RMSE3=sqrt(mean((obs-pr)^2)))
## [1] 729.864
(MAE3=mean(abs(obs-pr)))
## [1] 681.6363
(RMSPE3=sqrt(mean(((obs-pr)/obs)^2)))
## [1] 0.03018659
(MAPE3=mean(abs(obs-pr)/obs))
## [1] 0.02830442
(CI3=mean(ul-ll))
## [1] 2541.563
Errores más grandes que las anteriores predicciones
Resultados
res=data.frame(
par=c(length(coef(mod1)),length(coef(modEC)),length(coef(modEClin))+nrow(mod.atip$atip)),
sigma2=c(mod1$sigma2,modEC$sigma2,modEClin$sigma2),
AIC=c(AIC1,AIC2,AIC3),
BIC=c(BIC1,BIC2,BIC3),
RMSE=c(RMSE1,RMSE2,RMSE3),
MAE=c(MAE1,MAE2,MAE3),
RMSPE=c(RMSPE1,RMSPE2,RMSPE3),
MAPE=c(MAPE1,MAPE2,MAPE3),
CIml=c(CI1,CI2,CI3)
)
row.names(res)=c("ARIMA","ARIMA+EC","ARIMA+EC+OutTreat")
res
## par sigma2 AIC BIC RMSE MAE
## ARIMA 3 0.0009084550 -1370.715 -1355.471 517.8877 434.0531
## ARIMA+EC 5 0.0008518283 -1388.719 -1365.852 585.6177 461.5682
## ARIMA+EC+OutTreat 25 0.0004872871 -1540.034 -1444.068 729.8640 681.6363
## RMSPE MAPE CIml
## ARIMA 0.02119971 0.01784540 3441.743
## ARIMA+EC 0.02357383 0.01882707 3384.447
## ARIMA+EC+OutTreat 0.03018659 0.02830442 2541.563